// Test code for one sector 
// code for multi sector, multi country, optimal tax, tariff, and industrial policy
// Markov case
// Three countries, home country is 1, foreign + ROW
// transition from high trade cost to low trade cost during transition path
//Note is at endT_Msector_Mcountry_withTariff_externality.pdf, page 22 --

%========================================================================%
% Preamble
%========================================================================%
@#define countries = [1,2]
@#define countries_no1 = [2]
@#define sectors = [1]

var 
@#for co in countries
	x@{co},w@{co},riskfree@{co},
	@#if co!=1
		mu@{co},P@{co},
	@#endif
	@#for so in sectors
		T@{co}@{so},Lp@{co}@{so},Lr@{co}@{so},M@{co}@{so},GT@{co}@{so},rr@{co}@{so},
		@#if co!=1
			gammar@{co}@{so},taux@{co}@{so},taf@{co}@{so},gammaT@{co}@{so},
		@#endif
		@#for co2 in countries
			pi@{co}@{co2}@{so},
			@#for so2 in sectors
                @#if co!=1
				dGT@{co}@{so}@{co2}@{so2},
                @#endif
			@#endfor
		@#endfor
	@#endfor
@#endfor
@#for so in sectors
    wedge@{so},
@#endfor
uc,sa11
;

% exogenous variables
varexo ea11;

parameters beta,theta,gL, sigma, rho,dbar,delta,gammaL21,gammaL22,
@#for co in countries
    L@{co}
    @#for so in sectors
        alpha@{co}@{so}
        @#if co!=1
            @#for co2 in countries
                @#for so2 in sectors
                    gcoef0_@{co}@{so}@{co2}@{so2}
                @#endfor
            @#endfor
        @#endif
    @#endfor
@#endfor
@#for so in sectors
    betas@{so},
@#endfor
rhoa,std_ea,ncountries,nsectors,rbest,P1,rhod,std_ed,dcost11,dcost22,dcost12,dcost21,sa21;



%========================================================================%
% Parameter values
%========================================================================%

load setpara_stoch.mat

ncountries     =  mpara(1); % number of countries
nsectors       =  mpara(2);
beta   = mpara(3);     % discount factor
theta  = mpara(4);
sigma = mpara(5);
gL      = mpara(6);
dbar   = mpara(7);
delta = mpara(8);
rho   = mpara(9);
rbest = delta/(theta*(1-beta*(1-delta))+delta);

P1       = 0;
gammaL21 = 0;
gammaL22 = 0;
sa21=0;


ipara = 9;
  
@#for co in countries
	ipara = ipara+1;
	L@{co}=mpara(ipara);
@#endfor

@#for co in countries
	@#for so in sectors
		ipara = ipara+1;
		alpha@{co}@{so} = mpara(ipara);
	@#endfor
@#endfor 	

@#for so in sectors
	ipara = ipara+1;
	betas@{so} = mpara(ipara);
@#endfor

ipara = ipara+1;
rhoa   = mpara(ipara);
ipara = ipara+1;
std_ea = mpara(ipara);

ipara = ipara+1;
rhod   = mpara(ipara);
ipara = ipara+1;
std_ed = mpara(ipara);

@#for co in countries_no1
    @#for so in sectors
        @#for co2 in countries
            @#for so2 in sectors
    	        ipara = ipara+1;
                gcoef0_@{co}@{so}@{co2}@{so2}= mpara(ipara);   
            @#endfor
        @#endfor
    @#endfor
@#endfor


@#for co in countries
    @#for co2 in countries
		@#if (co==co2)
			dcost@{co}@{co2}=1.0;
		@#else
			dcost@{co}@{co2}=dbar;
		@#endif
	@#endfor
@#endfor


%========================================================================%
model;
%========================================================================%
// define uc

uc = exp(x1)^(-sigma);


[name='dG/dT guess']
@#for co in countries_no1 //destination
   @#for so in sectors
      @#for co2 in countries 
        @#for so2 in sectors   
            dGT@{co}@{so}@{co2}@{so2}= gcoef0_@{co}@{so}@{co2}@{so2}*(exp(x@{co}(+1))^(-sigma)*exp(P@{co}(+1))^(sigma-1)*exp(w@{co}(+1))/(exp(sa@{co}@{so}(+1))*alpha@{co}@{so}))/exp(T@{co2}@{so2});
         @#endfor
      @#endfor
   @#endfor
@#endfor

[name='Defintion: GT']
@#for co in countries
	@#for so in sectors
    	exp(GT@{co}@{so})= exp(x@{co})^(-sigma)*exp(P@{co})^(sigma-1)*exp(w@{co})/(exp(sa@{co}@{so})*alpha@{co}@{so});
    @#endfor
@#endfor


[name='Defintion: M']
@#for co in countries //destination
	@#for so in sectors
		M@{co}@{so} = exp(x@{co}(+1))^(-sigma)*exp(P@{co}(+1))^(sigma-1)*exp(w@{co}(+1))/(exp(sa@{co}@{so}(+1))*alpha@{co}@{so})*exp(x@{co})^sigma*exp(P@{co})^(1-sigma);
	@#endfor
@#endfor


[name='Pi matrix']
@#for co in countries //destination
    @#for co2 in countries 
        @#for so in sectors
            @#if (co2==1)
                @#if (co==1)
                    pi@{co}@{co2}@{so}=exp(T@{co2}@{so})*(exp(w@{co2}))^(-theta)/
                    	(exp(T@{co2}@{so})*(exp(w@{co2}))^(-theta)
                    		@#for co3 in countries_no1 //n
                    			+exp(T@{co3}@{so})*(exp(w@{co3})*(1+taf@{co3}@{so})*dcost@{co}@{co3})^(-theta)
                    		@#endfor
                    	)
                    ;
                @#else
                    pi@{co}@{co2}@{so}=exp(T@{co2}@{so})*(exp(w@{co2})*(1+taux@{co}@{so})*dcost@{co}@{co2})^(-theta)/
                        (exp(T1@{so})*(exp(w1)*(1+taux@{co}@{so})*dcost@{co}@{co2})^(-theta)
                            @#for co3 in countries_no1
                            	+exp(T@{co3}@{so})*(exp(w@{co3})*dcost@{co}@{co3})^(-theta)
                            @#endfor
                        )
                    ;
                @#endif
            @#else
                 @#if (co==1)
                     pi@{co}@{co2}@{so}=exp(T@{co2}@{so})*(exp(w@{co2})*(1+taf@{co2}@{so})*dcost@{co}@{co2})^(-theta)/
                           (exp(T1@{so})*(exp(w1))^(-theta)
                                @#for co3 in countries_no1
                                   +exp(T@{co3}@{so})*(exp(w@{co3})*(1+taf@{co3}@{so})*dcost@{co}@{co3})^(-theta)
                                @#endfor
                            )
                        ;
                    @#else
                        pi@{co}@{co2}@{so}=exp(T@{co2}@{so})*(exp(w@{co2})*dcost@{co}@{co2})^(-theta)/
                            (exp(T1@{so})*(exp(w1)*(1+taux@{co}@{so})*dcost@{co}1)^(-theta)
                                @#for co3 in countries_no1
                                    +exp(T@{co3}@{so})*(exp(w@{co3})*dcost@{co}@{co3})^(-theta)
                                @#endfor
                            )
                        ;
                    @#endif
            @#endif 
       @#endfor //so
    @#endfor //co2
@#endfor //co, destination
        

[name='T evolution']
@#for co in countries
    @#for so in sectors
        exp(T@{co}@{so})- (1-delta)*exp(T@{co}@{so}(-1)) = (exp(sa@{co}@{so})*alpha@{co}@{so})*exp(Lr@{co}@{so});
    @#endfor
@#endfor


[name='Country 1 price=1']
@#for  so in sectors
    -betas@{so}/theta*
    log(exp(T1@{so})*(exp(w1))^(-theta)+
        @#for co2 in countries_no1
           +exp(T@{co2}@{so})*(exp(w@{co2})*(1+taf@{co2}@{so})*dcost@{1}@{co2})^(-theta)
        @#endfor
     )
@#endfor
=0;

[name='Country n price']
@#for co in countries_no1
    @#for so in sectors
        -betas@{so}/theta*
        log(
            exp(T1@{so})*(exp(w1)*(1+taux@{co}@{so})*dbar)^(-theta)
            @#for co2 in countries_no1
               +exp(T@{co2}@{so})*(exp(w@{co2})*dcost@{co}@{co2})^(-theta)
            @#endfor
           )
    @#endfor
    =P@{co};
@#endfor


[name='market clearing C1']
@#for so in sectors
    (1+theta)/theta*exp(w1)*exp(Lp1@{so})=betas@{so}*
            (
                pi11@{so}*exp(x1)
                @#for co in countries_no1
                    +pi@{co}1@{so}*exp(x@{co})/(1+taux@{co}@{so})
                @#endfor
            )
            ;
@#endfor


//[name='labor market clearings, N']

@#for co2 in countries
    @#for so in sectors
        +exp(Lr@{co2}@{so})+exp(Lp@{co2}@{so})
    @#endfor
    =L@{co2};
@#endfor

[name='expenditure country 1: x1']
exp(x1)= (1+theta)/theta*exp(w1)*
        (
        @#for so in sectors
        	+ exp(Lp@{1}@{so})
        @#endfor
        )
        @#for co in countries_no1
        	@#for so in sectors
        		+ betas@{so}*taux@{co}@{so}/(1+taux@{co}@{so})*pi@{co}@{1}@{so}*exp(x@{co})
        	@#endfor
        @#endfor
        @#for co in countries_no1
        	@#for so in sectors
            	+betas@{so}*taf@{co}@{so}/(1+taf@{co}@{so})*pi@{1}@{co}@{so}*exp(x1)
            @#endfor
        @#endfor
    ;

[name='expenditure country n: xn']
@#for co in countries_no1
    exp(x@{co}) = (1+theta)/theta*exp(w@{co})*
    	(
    		@#for so in sectors
    			+ exp(Lp@{co}@{so})
    		@#endfor
    	)
    	;
@#endfor


[name='FOC on TNj, Ns*(N-1)']

@#for co in countries_no1 //n
@#for so in sectors
    -uc*exp(x1)*betas@{so}*pi@{1}@{co}@{so}
    +theta*uc*(
         @#for co2 in countries_no1 //m
            +betas@{so}*taux@{co2}@{so}/(1+taux@{co2}@{so})*pi@{co2}1@{so}*pi@{co2}@{co}@{so}*exp(x@{co2})
          @#endfor
        )
    -theta*beta*exp(T@{co}@{so})*gammaT@{co}@{so}(+1)*(1-delta)
   = gammaL@{co}@{so}*theta*betas@{so}*(
         @#for co3 in countries_no1 //m
             +pi@{co3}@{co}@{so}*exp(x@{co3})
          @#endfor
        )
   +(
        @#for co2 in countries_no1 //i
          -gammaL@{co2}@{so}*theta*betas@{so}*(
            @#for co3 in countries_no1 //m
               +pi@{co3}@{co2}@{so}*pi@{co3}@{co}@{so}*exp(x@{co3})
           @#endfor //co3
         )
       @#endfor //co2
     )
    - gammar@{co}@{so}*exp(w@{co})*exp(Lp@{co}@{so})/exp(T@{co}@{so})
    - (
        @#for co2 in countries_no1 //i
            @#for so2 in sectors //k
             + gammar@{co2}@{so2}*beta*(1-delta)*(1-sigma)*M@{co2}@{so2}*betas@{so}*pi@{co2}@{co}@{so}
            @#endfor
        @#endfor
      )
    +(
        @#for co2 in countries_no1 //i
           @#for so2 in sectors //k
             + theta*gammar@{co2}@{so2}*beta*(1-delta)*dGT@{co2}@{so2}@{co}@{so}*exp(x@{co2})^(sigma)*exp(P@{co2})^(1-sigma)*exp(T@{co}@{so})
           @#endfor
       @#endfor
     )
    - theta*exp(T@{co}@{so})*gammaT@{co}@{so}
    ;
@#endfor //so
@#endfor

[name='FOC on Lpj, Ns*(N-1)']

@#for co in countries_no1 //n
    @#for so in sectors //j
      -uc*(
          @#for so2 in sectors //k
            +betas@{so2}*taux@{co}@{so2}/(1+taux@{co}@{so2})*pi@{co}1@{so2}
          @#endfor
         )
      = -gammaL@{co}@{so}
        +(
            @#for co2 in countries_no1 //i
               @#for so2 in sectors //k
                 +gammaL@{co2}@{so2}*betas@{so2}*pi@{co}@{co2}@{so2}
               @#endfor
            @#endfor
         )
        + gammar@{co}@{so}/exp(T@{co}@{so})/(1+theta)
        + (          
           @#for so2 in sectors //k
             +gammar@{co}@{so2}*beta*(1-delta)*sigma*M@{co}@{so2}/exp(x@{co})
           @#endfor
         )
        - theta/(1+theta)*mu@{co}/exp(w@{co})
      ;  
    @#endfor //so
@#endfor


[name = 'FOC on Lr with Ns(N-1)']

@#for co in countries_no1
    @#for so in sectors
        gammaT@{co}@{so}*(exp(sa@{co}@{so})*alpha@{co}@{so}) = mu@{co};
    @#endfor
@#endfor

[name='FOC for wn, n>2']

@#for co in countries
@#if (co>2)
    uc*exp(x1)*(
        @#for so in sectors 
            +betas@{so}*pi1@{co}@{so}
        @#endfor
       )
    -uc*(
        @#for so in sectors
        	+betas@{so}*taux@{co}@{so}/(1+taux@{co}@{so})*pi@{co}1@{so}*exp(x@{co})
        	+betas@{so}*(
        		@#for co2 in countries_no1
        			+taux@{co2}@{so}/(1+taux@{co2}@{so})*theta*pi@{co2}1@{so}*pi@{co2}@{co}@{so}*exp(x@{co2})
        		@#endfor
        	)
        @#endfor
     )
    =(
    	@#for so in sectors
    		-gammaL@{co}@{so}*(1+theta)/theta*exp(w@{co})*exp(Lp@{co}@{so})
    		-gammaL@{co}@{so}*theta*betas@{so}*(
    			@#for co2 in countries_no1
    				+pi@{co2}@{co}@{so}*exp(x@{co2})
    			@#endfor
    		)
    	@#endfor
    )
	+(
		@#for co2 in countries_no1
		@#for so in sectors
			+gammaL@{co2}@{so}*betas@{so}*(
				@#for co3 in countries_no1 
					+theta*pi@{co3}@{co2}@{so}*pi@{co3}@{co}@{so}*exp(x@{co3})
				@#endfor
				+ pi@{co}@{co2}@{so}*exp(x@{co})
			)
		@#endfor
		@#endfor
	)
	+(
		@#for so in sectors 
			+gammar@{co}@{so}*beta*(1-delta)*(sigma-1)*M@{co}@{so}
		@#endfor
	)
	+(
		@#for co2 in countries_no1 
		@#for so in sectors 
			+gammar@{co2}@{so}*beta*(1-delta)*(1-sigma)*M@{co2}@{so}*(
				@#for so2 in sectors 
					+betas@{so2}*pi@{co2}@{co}@{so2}
				@#endfor
			)
		@#endfor
		@#endfor
	)
	;
@#endif
@#endfor

//[name='innovation choice for home country']

@#for so in sectors
    exp(w1)/(exp(sa1@{so})*alpha1@{so})
        =  1/theta*exp(w1)*exp(Lp1@{so})/exp(T1@{so})+beta*(1-delta)*M1@{so}
          + theta/(1+theta)/uc*
              (
                @#for co in countries_no1 //n
                    @#for so2 in sectors
                        +gammar@{co}@{so2}*beta*(1-delta)*dGT@{co}@{so2}1@{so}*exp(x@{co})^sigma*exp(P@{co})^(1-sigma)
                     @#endfor
                @#endfor
               )
     ;      
@#endfor

//[name='wedge']

@#for so in sectors
    wedge@{so}
        = theta/(1+theta)/uc*
              (
                @#for co in countries_no1 //n
                    @#for so2 in sectors
                        +gammar@{co}@{so2}*beta*(1-delta)*dGT@{co}@{so2}1@{so}*exp(x@{co})^sigma*exp(P@{co})^(1-sigma)
                     @#endfor
                @#endfor
               )/(1/theta*exp(w1)*exp(Lp1@{so})/exp(T1@{so}))
     ;      
@#endfor

[name=' optimal export tax , Ns*(N-1)']
@#for co in countries_no1
    @#for so in sectors
        1+taux@{co}@{so} = uc*(1+theta*(1-pi@{co}1@{so}))/
        (uc*theta*(1-pi@{co}1@{so})
        	-(
        		@#for co2 in countries_no1 //i
        			+theta*gammaL@{co2}@{so}*pi@{co}@{co2}@{so}
        		@#endfor
        	)
        	+ (
        		@#for so2 in sectors //k
        			+gammar@{co}@{so2}*beta*(1-delta)*(sigma-1)*M@{co}@{so2}/exp(x@{co})
        		@#endfor
        	)
        );
    @#endfor
@#endfor

[name='optimal tariff, Ns*(N-1)']
@#for co in countries_no1
    @#for so in sectors 
        1+taf@{co}@{so} = (uc-gammaL@{co}@{so})/uc;
    @#endfor
@#endfor

[name=' innovation choice for non-home countries']
@#for co in countries_no1
    @#for so in sectors
        exp(w@{co})/(exp(sa@{co}@{so})*alpha@{co}@{so})=1/theta*exp(w@{co})*exp(Lp@{co}@{so})/exp(T@{co}@{so})+
                beta*(1-delta)*M@{co}@{so}
        ;
    @#endfor
@#endfor

@#for co in countries
    @#for so in sectors
    	exp(rr@{co}@{so}) = exp(Lr@{co}@{so})/(exp(Lr@{co}@{so})+exp(Lp@{co}@{so}));
	@#endfor
@#endfor

[name='defintion of risk free rate']
@#for co in countries    
    exp(riskfree@{co}) = exp(x@{co})^(-sigma)*exp(P@{co})^(1-sigma)/(beta*exp(x@{co}(+1))^(-sigma)*exp(P@{co}(+1))^(1-sigma));
@#endfor

[name='shocks']
sa11 = rhoa*sa11(-1)+ std_ea*ea11;



end;



//****************************************************************************

steady;
check;

shocks;
//@#for co in countries
//    @#for so in sectors
//    	var ea@{co}@{so} = 1;
//	@#endfor
//@#endfor
var ea11=1;
//var ed = 1;
//var ea12; stderr std_ea;
//var ea21; stderr std_ea;
//var ea22; stderr std_ea;
end;

stoch_simul(order=1,irf=200,nograph,noprint) sa11,T11,T21,taux21,rr11,rr21,x1,x2,riskfree1,riskfree2,wedge1,w1,w2,P2;

//sa1,T1,T2,taux2,r1,r2,x1,x2,riskfree1,riskfree2;
//T11,T21,taux21,taf21,sa11,Lr11,Lr21,rr11,rr21,wedge1,w1,w2,x1,riskfree1,riskfree2,x2;